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ABSTRACT 

Context. Supermassive black holes are probably present in the centre of the majority of the galaxies. There is a consensus 
that these exotic objects are formed by the growth of seeds either by accreting mass from a circumnuclear disk and/or 
by coalescences during merger episodes. 

Aims. The mass fraction of the disk captured by the central object and the related timescale are still open questions, 
as well as how these quantities depend on parameters like the initial mass of the disk or the seed or on the angular 
momentum transport mechanism. This paper is addressed to these particular aspects of the accretion disk evolution 
and of the growth of seeds. 

Methods. The time-dependent hydrodynamic equations were solved numerically for an axi-symmetric disk in which the 
gravitational potential includes contributions both from the central object and from the disk itself. The numerical code 
is based on a Eulerian formalism, using a finite difference method of second-order, according to the Van Leer upwind 
algorithm on a staggered mesh. 

Results. The present simulations indicate that seeds capture about a half of the initial disk mass, a result weakly 
dependent on model parameters. The timescales required for accreting 50% of the disk mass are in the range 130- 
540 Myr, depending on the adopted parameters. These timescales permit to explain the presence of bright quasars at 
2 ~ 6.5. Moreover, at the end of the disk evolution, a "torus-like" geometry develops, offering a natural explanation 
for the presence of these structures in the central regions of AGNs, representing an additional support to the unified 
model. 
Conclusions. 
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1. Introduction 

At the present time, the general belief is that supermas- 
sive black holes (SMBHs) located at the center of galax- 
ies have been formed by the growth of primordial "seeds 
cither by matter accretion or coalescences during merger 
episodes. This picture is consistent with the fact that the 
present black hole (BH) mass density compares quite well 
with the mass density derived from the bolometric lumi- 
nosity function of quasars, under the assumption that the 
accretion process itself is the source of the radiated energy 
(Soltan 1982; Small & Blandford 1992; Hopkins, Richards 
& Hernquist 2007). The accreted matter is essentially bary- 
onic in origin with an eventual contribution of the exotic 
dark component not exceeding 10% of the total accreted 
matter (Peirani & de Freitas Pacheco 2008). 

Nevertheless, the formation of massive BHs by the di- 
rect collapse of central regions of proto-galaxies is an alter- 
native possibility considered by different authors. Shapiro 
(2004) considered the core collapse of relativistic star clus- 
ters formed in starbursts, which may have occurred in the 
early evolution of galaxies. The collapse of a collisionlcss 
system, after having experienced a phase of strong oscil- 
lations followed by violent relaxation, may undergo a slow 
gravothermal instability, if the phase space incomprcssibil- 
ity condition is violated (Levin, Pakter & Rizzato 2008), 



leading eventually to the formation of a BH. Despite the 
non inclusion of the star formation process and its asso- 
ciated feedback, Wise, Turk & Abel (2009) found from 
hydrodynamical simulations of <~ 10 s M proto-galaxies, 
that massive central objects with masses around 10 5 M Q 
could formed inside regions of one parsec radius. The most 
favourable physical conditions for the gas in these proto- 
halos to form those massive central objects are when the 
virial temperature is about 15000 K and the circular veloc- 
ity is around 20 km/s (Reagan & Haehnelt 2009). 

If the aforementioned investigations open the possibil- 
ity of forming directly massive BH via the gravitational 
collapse, other studies lead to the formation of interme- 
diate mass BH, which could be the required seeds. Black 
holes in the mass range (10 3 — 10 4 M Q ) could be formed 
in the collapse of primordial gas clouds (Haehnelt & Rees 
1993; Eisenstein & Loeb 1995; Koushiappas, Bullock & 
Dekel 2004). Another appealing possibility, which will be 
considered in the present investigation, is the formation of 
BH seeds with masses in the range 10 2 — 10 4 Mq, result- 
ing from the evolution of primordial massive stars (Bromm, 
Coppi & Larson 1999; Abel, Bryan & Norman 2000; Heger 
et al. 2003; Yoshida et al. 2006). These stars would have 
been formed at z <~ 15 — 20 in the high density peaks of 
the primordial fluctuation spectrum (Gao et al. 2007) and 
their high masses would be the consequence of very inef- 
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ficient gas cooling at zero metallicity. Besides forming BH 
seeds, these massive stars could also have contributed to 
the reionization of the universe. 

Different investigations based on cosmological simula- 
tions have been performed in these past years, aiming to un- 
derstand not only the growth process of seeds, but also the 
consequences of the feedback resulting from the accretion 
process itself on the environment as well as the observed 
correlations between the black hole mass with properties of 
the host galaxy (DiMatteo et al. 2003; Pelupessy, DiMatteo 
Ciardi 2007; Springel, DiMatteo & Hernquist 2005; Filloux 
et al. 2009). In all these simulations, the physical descrip- 
tion of the accretion process requires a substantial simpli- 
fication, since it occurs on unresolved spatial scales. If the 
BH is at rest with respect to the gas, the inflow is probably 
spherical, almost adiabatic and the accretion rate is esti- 
mated by the Hoyle-Lyttlcton-Bondi (HLB) formula (Hoylc 
& Lyttleton 1939; Bondi & Hoyle 1944), whose rate depends 
on parameters like the gas density and the sound velocity 
evaluated far away from the BH, on scales supposed to be 
resolved by simulations. The assumption that BHs could be 
at rest and that the accretion process is spherically sym- 
metric is probably unrealistic but frequently used in cos- 
mological simulations due to its simplicity. However, the 
radiation emitted during the inflow, in particular near the 
BH horizon, affects the surrounding gas, reducing consider- 
ably the accretion rate and rendering inefficient the growth 
of seeds by such a mechanism (Milosavljevic et al. 2009). 
The situation is rather different if the BH is moving with 
respect to the gas. After passing the BH, a conically shaped 
shock is produced in the flow, in which the gas loses the mo- 
mentum component perpendicular to the shock front. After 
compression in the shock, gas particles within a certain im- 
pact parameter will fall into the black hole. Particles hav- 
ing angular momentum exceeding 2r s c, with r g being the 
gravitational radius, will form a disk and only after viscous 
stresses have transported away the excess of angular mo- 
mentum will the gas cross the BH horizon. Clearly, all these 
processes are not adequately described in any of the afore- 
mentioned simulations and one may wonder if the accre- 
tion or the resulting luminosity are estimated satisfactorily 
when the Hoyle-Lyttlcton-Bondi approach is adopted. 

Numerical simulations indicate that after the merger of 
two galaxies, a considerable amount of gas is settled into 
the central region of the resulting object. The gas loses 
angular momentum in a timescale comparable to the dy- 
namical timescale, forming a circumnuclear self-gravitating 
disk (Mihos & Hernquist 1996; Barnes 2002), having masses 
in the range 10 6 — 10 9 M Q and dimensions of about 
100 — 500 pc. Such a disk will probably be able to feed 
the central BH, increasing its mass. In fact, this scenario 
was considered by Filloux et al. (2009) in their simulations. 
These authors adopted an accretion rate derived from a 
steady disk model, in which the angular momentum trans- 
port is provided by turbulent viscosity. Since in reality the 
disk is expected to be non-steady, they assumed that dur- 
ing a time step the properties of the disk do not change 
appreciable but they are updated continuously at each new 
time step. Despite this simple description of the accretion 
process, they were able to reproduce quite satisfactory the 
main properties of BHs and their host galaxies. More re- 
cently, Power, Nayakshin & King (2010) introduced a new 
method to model the physics of an accreting disk around a 
BH, which has some common points with the approach by 



Filloux et al. (2009). Based on a series of numerical exper- 
iments, they have shown that the disk accretion mode is 
more physically consistent than the simple use of the HLB 
formula and more efficient to feed the seed, in agreement 
with the conclusions by Filloux et al. (2009). In these last 
two approaches, the accretion disk is represented in sim- 
ulations by a "single" (eventually two in rare cases) SPH 
particle and the adopted prescription for the accretion rate 
defines the mass fraction transferred to the central BH. In 
the real world, the structure of the disk changes contin- 
uously for different reasons: mass is transferred either to 
central BH or to the outskirts of the disk, stars form con- 
suming part of gas available to feed the BH and infall of 
matter is a source of fresh gas replenishing the disk. 

The understanding of the aforementioned aspects re- 
quire a close inspection of the accretion process, in order to 
improve the modelization adopted in cosmological simula- 
tions. In the present paper we report results on the numer- 
ical solution of the hydrodynamic equations describing the 
evolution of a circumnuclear disk around a black hole, aim- 
ing to answer some specific questions that could be useful 
to improve the modelling of such a process in cosmological 
simulations. Moreover, SLOAN quasars associated to very 
massive BHs (Fan et al. 2001; Willot, McLure & Jarvis 
2003) are already observed at redshifts z ~ 6, when the 
universe was only 0.96 Gyr old. These observation imply 
that the growth process of seeds should occur in shorter 
timescales, in regions where the accretion mechanism must 
be very efficient. The present investigation will be mainly 
focused in the following questions: 1) what is the fraction of 
the disk mass accreted by the central BH? 2) Does such a 
fraction depend on the initial mass either of the BH seed or 
the disk? 3) What are the typical accretion timescales and 
how they depend on the disk parameters? This paper is or- 
ganized as follows: in Section 2 the accretion disk model is 
discussed and the main equations are introduced; in Section 
3 the numerical methods are presented; in Section 4 the 
main results are given and, finally, in Section 5 the conclu- 
sions are summarized. 

2. The accretion disk model 

In this work, we will study the evolution of nuclear disks 
formed possibly after a merger event or by gas infall, hav- 
ing masses in the range 10 7 - 10 8 M & and with typical 
dimensions of the order of 50 pc. These disks are initially 
self-gravitating and, as the central BH grows, it dominates 
the dynamics of the inner regions while the outer parts are 
still dominated by the disk self gravitation. As we shall see 
later, the external disk is constituted by neutral gas having 
temperatures of about 100 — 2000 K, conditions favourable 
to form stars. However, in the present work, the gas con- 
version into stars through the disk will not be considered, 
although a new version of our code, presently under devel- 
opment, will offer such a possibility. 

2.1. Angular momentum transport 

A major problem in the current understanding of accretion 
disks is the absence of a physical theory able to describe the 
viscosity of the gas in the presence of turbulent flows or in 
the presence of a magnetic field. The angular momentum 
transport in accretion disks is generally described by the 
formalism introduced almost forty years ago by Shakura & 
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Sunyaev (1973), in which the viscosity due to the subsonic 
turbulence is parametrized by the relation 



r\ = aHc s 



where a < 1 is a dimensionless coefficient, H is the vertical 
scale of the disk, supposed to be of the same order as the 
typical (isotropic) turbulence scale £ t and c s is the sound 
velocity. Accretion disk models based on the "a-viscosity" 
approach are able to explain successfully the observed prop- 
erties of dwarf-novae and X-ray binaries, but disks based on 
such a formulation are thermally unstable as demonstrated 
long time ago by Piran (1978). 

The possibility that turbulence could be generated by 
local gravitational instability in geometrically thin disks 
was considered by Paczynski (1978). Duschl and Britsch 
(2006) revisited this idea by analysing the gravitational 
instability of self-gravitating disks, concluding that such 
an instability could be able to develop turbulence in the 
flow and, consequently, to generate viscosity without requir- 
ing the contribution of magneto-hydrodynamic turbulence. 
Recent simulations of the gas inflow in the central regions 
of galaxies, induced by the gravitational potential either of 
the stellar nucleus or the SMBH, reveal the appearance of 
highly supersonic turbulence, with velocities of the order 
of the virial value (Regan & Haehenelt 2009; Levine et al. 
2008; Wise, Turk & Abel 2008). Amazingly, no fragmenta- 
tion is observed in such a gas despite of being isothermal 
and gravitationally unstable. In fact, Begelman & Shlosman 
(2009) argued that an efficient angular momentum transfer 
suppresses fragmentation, favouring the gas inflow. On the 
contrary, if the angular momentum transfer is inefficient, 
the turbulence decays and triggers global instabilities which 
regenerates a turbulent flow. Moreover, according to their 
analysis, fragmentation is suppressed whenever the gas tem- 
perature remains below the virial value. As we shall see 
later, during the early evolutionary phases of the disk, the 
accretion rate is quite important. Consequently, the heat 
generated by dissipation of turbulence increases the radi- 
ation pressure which inflates the inner regions, producing 
a "slim" disk. We would expect that a flow self-regulated 
by the aforementioned mechanism could be stable against 
fragmentation. If the flow is self-regulated, it must be char- 
acterized by a critical Reynolds number 1Z, determined by 
the viscosity below which the flow becomes unstable. In 
fact, de Freitas Pacheco & Steiner (1976) suggested that in- 
stead of the "a" parametrization, the effective (turbulent) 
viscosity r\ could be expressed in terms of such a critical 
Reynolds number by 



V = 



n 



(2) 



where r is the radial distance to the center of the disk (or 
to the central BH) and V$ — fir is the azimuthal velocity of 
the gas. It is worth mentioning that this formulation, which 
will be adopted here, is essentially the "/3-viscosity" model 
discussed by Duschl, Strittmatter & Biermann (1998) (see 
also Richard & Zahn 1999). An additional aspect justify- 
ing the adoption of this formulation is that accretion disks 
modeled by such a viscosity prescription are thermally sta- 
ble according to the analysis by Piran (1978). However, 
in the present work the disk stability was not investigated 
and we cannot exclude the possibility that fragmentation 
occurs. This aspect will be examined in a future paper. 



It can be verified trivially that in such a formulation 
the Mach number associated to the mean turbulent mo- 
tions is M t ~ V$/(c s yR,). This implies that turbulence is 
slightly supersonic in most regions of the disk but it can be 
supersonic at the inner regions, depending on the critical 
Reynolds number. 

2.2. The dynamical equations 

The hydrodynamic equations are written in cylindrical co- 
ordinates (r, <p, z) and, in order to simplify the solution of 
the Poisson equation, variables are integrated along the z- 
axis but the scale of height is calculated in a consistent way 
as we shall see below. The equation of the radial motion, 
neglecting gas pressure gradients, is 



dt r dr r dr 







(3) 



where V r is the radial velocity, is the azimuthal velocity 
and tp is the gravitational potential, including the contri- 
bution of the central BH and of the disk itself. The central 
BH is probably rotating but in the present investigation we 
discard this possibility and we use the approximate poten- 
tial of Paczynski- Wiita (Paczynski & Wiita 1980) to model 
the gravitational field of the BH. The contribution of the 
disk itself is obtained from the expression for the internal 
potential of an ellipsoid in which we have performed the 
limit z — > 0. Under these conditions the total gravitational 
force is 

dip _ GM BH {t) G f r E(a,t)ada 
dr 



(r - r g ) 2 



o Vr 



(4) 



where r g = 2GMbh /c 2 is the gravitational (or the event 
horizon) radius and E is the columnar mass density defined 
by 



E(r, t) = / p(r,z,t)dz 

J — oo 

The mass of the BH varies according to the equation 
dM BH {t) 



dt 



= 2irri so (t)Y,(ri so ,t)V r (ri so ,t) 



(5) 



(6) 



where the quantities on the right side are evaluated at the 
radius of the last stable circular orbit ri so . Notice that as the 
BH grows, the radius of the horizon and of the last stable 
circular orbit increase proportionally to the BH mass. 
The continuity equation is given by 

<9E 1 d(rHV r 
~8t * 







(7) 



r dr 

and the equation for the azimuthal motion, including the 
viscous forces responsible for the angular momentum trans- 
port is 



dV<p | v dV^ | VrVt | 1 d(rT r<j> ) 
dt r dr r rE dr 







(8) 



where the considered component of the stress tensor (inte- 
grated along the z-axis) is 



T r4> = ryEr 



dtt 

dr 



and, using eq. [2j the equation above can be recast as 



T, 



T(j> 



TZ dr 



(9) 



(10) 
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2.3. The scale of height 

Assuming that the disk is in hydrostatic equilibrium along 
the z-axis, we can write 



d(P g +P t + P r ) 

dz 



= ~P9z 



(11) 



where P g = pc 2 s , P t = p < v 2 > and P r = aT 4 /3 are 
respectively the gas, the turbulent and the radiation pres- 
sure. The vertical component of the gravitational acceler- 
ation has two contributions: one from the disk self-gravity 
and another from the tidal force due to the central black 
hole. Thus, 



o r^v z , GMbh 
g z = 27rGE— H = — z 

H r A 



(12) 



Defining an effective scale of height by the relation H = 
p/(\ dp/dz |) and using the equations above, one obtains 
after some algebra 



H = — 



A" 



(l-/3)-V(l-/3) 2 + Q 2 (l + e 2 ) 



(13) 



where e 2 =< v\ > /c~ and we have introduced respectively 
GM BH 



O 2 



and 



Q = 



VL K c s 



ttGE 

The parameter (3 defined as 



3ttGS 2 



measures the ratio between the radiation pressure and the 
disk self-gravity. Notice that from eq.[l3]some limiting cases 
can be derived. When the radiation force is unimportant 
(f3 << 1) and self-gravity dominates (Q << 1), one obtains 



H 



c 2 (l + £ 2 ) 
2ttGE 



(17) 



whereas when Q » 1, namely, the vertical acceleration is 
dominated by tidal forces we have 



H 



+ £ 2 ) 



(18) 



It is worth mentioning that when the disk becomes opti- 
cally thin, the contribution due to the radiation pressure is 
negligible and, in this case the scale of height is computed 
from eq. 13 with the condition /3=0. 



2.4. The temperature of the disk 

We assume that all energy is generated locally by vis- 
cous dissipation and radiated essentially along the z-axis 
although a small fraction could be advected inwards. Under 
these conditions, the local energy balance is given by 



; dis 



idv 



(19) 



The energy rate per unit of area dissipated by viscous forces 
is 

(20) 



Qau = %^r 2 ( 5^0 



while the advected energy flux is (Abramowicz et al. 1986) 

(21) 



Qadv = f(r) ( — ) (J, I,, 



where f(r) depends on the angular momentum distribution 
throughout the disk but, in general, it is of the order of the 
unity (Abramowicz et al 1995), except near the star-disk 
transition layer, where it could attain quite large values. 
In our case, since the central object is a black hole, such a 
layer doesn't exist and the f(r) factor depends essentially 
on the difference between the angular momentum at the 
considered position and that at the last stable orbit. Here 
we assume simply that f(r) = 1 everywhere in the disk. 

In the regions where the disk is partially or completely 
ionized, the opacity is due to the Thomson scattering, 
bound-free and free-free transitions. In this case, using 
the solution of the transfer equation within the Eddington 
approximation derived by de Freitas Pacheco, Steiner & 
Daminelli (1977), the local emergent flux is 



(14) $„(T e (r)) 



2nB v (T e 



1 



2A 



0,v 



(22) 



where to t „ = ySry (jf + t s ) is the effective optical depth, 
with Tf and t s being respectively the absorption and scat- 
tering optical depths. The other parameter is defined by the 
relation A 2 = Tf /(Tf + t s ). T e is the effective temperature. 
Thus, 



(23) 



(15) 



(16) Qrad = 2 / $„{T e {r))dv 



where the factor 2 takes into account both disk surfaces. 

Another factor affecting the evaluation of the disk tem- 
perature is the photon trapping. At very high optical 
depths, the photon diffusion timescale along the z-axis 
td = 3(H/c)to can be larger than the accretion timescale 
tac = f/V r . Here io is the frequency averaged effective opti- 
cal depth as defined above. When this occurs, photons are 
not able to reach the surface, being convected inwards and 
swallowed by the central BH. In this case, besides the ad- 
vection correction, the dissipated energy should be reduced 
by an extra factor equal to the ratio t ac /td- 

The local effective temperature is derived from the nu- 
merical solution of the above equations. In the outer re- 
gions, the disk is neutral and optically thin. In this case, 
the temperature was estimated simply from the balance be- 
tween the heating and the cooling rates. The cooling is due 
essentially to molecular hydrogen in the case of a disk con- 
stituted essentially by primordial gas. When trace elements 
are present, an additional cooling mechanism should be 
considered, which is due to the excitation of fine-structure 
levels of C + , Si + , O and Fe + . The chemical reactions lead- 
ing to the formation of the Hi molecule requires a residual 
electron density, which we will assume to be provided by 
the ionizing radiation background. Under this conditions, 
we have adopted the cooling functions for molecular hy- 
drogen given by Galli & Palla (1998) and those for atomic 
infrared lines by de Freitas Pacheco (1969). 
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3. The numerical method 

In order to solve numerically the hydrodynamic equations 
describing the structure and the evolution of the disk, we 
have modified the public code FARGO (Masset 2000) , origi- 
nally developed for studies of the interaction between plan- 
ets and circumstellar disks. In fact, the structure of this 
code is quite similar to ZEUS (Stone & Norman 1992). The 
code is based on an Eulerian formalism, using a finite dif- 
ference method of second-order according to the Van Leer 
upwind algorithm on a staggered- mesh (Van Leer 1977), 
which means that scalar quantities like surface density, scale 
of height, etc., are defined at the center of the cell whereas 
vector quantities like velocity, fluxes in general, etc., are de- 
fined at the interface between cells. In our particular case, 
the grid implementation was done in a logarithmic scale, 
since the radial scale may vary between nine up to eleven 
orders of magnitude in order to cover all the extension of 
the disk from the last stable orbit until the external radius. 
Such a logarithmic grid was covered by 1024 ring sectors. 

The time-step is controlled by the Courant-Friedrich- 
Levy (CFL) condition, which states that the information 
cannot sweep a distance larger than the size of a cell 
(Ri — Ri-i) over one time-step. This condition introduces 
five constraints, corresponding respectively to five time-step 
limits 8t\, St 2 , St 3 , Sti and St$. The first one is defined by 



Sti = min(Sr,rS<f))/c s 



(24) 



where Sr and rS<j) correspond to the radial and azimuthal 
mesh size. This condition implies that no wave propagat- 
ing with velocity c s should cross the cell in one time-step. 
The following two conditions constrain the motion of a test 
particle that should not travel a distance greater than Sr 
in the radial direction and rScf) in the azimuthal direction. 
These conditions are expressed as 



St 2 = Sr/V r 
and 

St 3 = rS<j>/V^ 



(25) 



(26) 



The fourth condition is related to the artificial viscosity, 
introduced to avoid large discontinuities due to the eventual 
presence of shocks (Richtmyer & Morton 1957) and is given 
by 



5t± = min 



Sr 



(ir i ic*dv„ 



(27) 



where SV r = d r V r Sr and 5V<f, = d^V^Scf). Notice that for ax- 
ial symmetry this last term is zero. C v is the Von Neumann- 
Richtmyer viscosity constant, usually taken in the range (1, 
2). Here, we have assumed C v — 1.41 in all runs. The last 
condition assures that two neighboring rings are not dis- 
connected after a time-step and it can be expressed as 



SU = 



2tt 



[Q(r) - Q(r + Sr)] 



(28) 



Finally, the hydrodynamic time-step is computed from the 
relation 



At = Ccfl y~min 



1 



St^ 



St 2 2 



-;St 5 



Stl 



■st: 



(29) 



where Ccfl is the so-called CFL number, usually taken in 
the range (0, 1) and assumed here to be equal 0.5 in all 
runs. 

In order to understand the different steps of the algo- 
rithm, we recall that the hydrodynamic equations can be 
written in a general form, i.e., 



d t X + VAX = S 



(30) 



where X is a generic variable (density, mass flux or veloc- 
ity) and S is a generic source term. Our algorithm pro- 
ceeds in the following way. After evaluation of the gravita- 
tional potential, the velocity field is updated in each cell, 
using the source terms and a first order integrator such as: 
X(t+dt)=X(t)+Sdt. Then, the considered variable is ad- 
vected from each cell to its neighbour with a second order 
upwind interpolation. The quantities S, T,V r and rSV^ are 
updated by taking into account their associated fluxes. 

The boundary conditions arc simple equations which as- 
sign values to the dependent variables in the ghost zones 
(inner and outer rings of the grid) corresponding to values 
in the adjacent ring. The fundamental rule applied in the 
ghost zones is that matter are not allowed to enter in the 
grid zone. Only outflows arc permitted. Thus, in the inner 
ghost zone (ring zero) , the density is that of the neighbour- 
ing ring (ring one). The radial velocity of the ring one is set 
equal to that of ring 2 if the later is negative (outflow) and 
equal to zero otherwise. A symmetric algorithm is used for 
the outer ghost zone. 

For the initial configuration, we assume that the colum- 
nar density has a profile given by 



£(r,t = 0) = £ - 
r 



(31) 



corresponding to an initial total disk mass Md = 2irYiQroRd, 
where Rd is the external radius of the disk. Supposing that 
the disk is initially in equilibrium, the considered density 
distribution implies that the initial profile of the azimuthal 
velocity is 



(32) 



while the initial profile of the radial velocity is simply 
V r (r,t = 0) =0. 

We have adopted the parsec as the unit of length, the 
initial BH mass as the unit of mass and the unit of time is 



GM BH (0) 



= 1.494xl0 6 ^^-J 



3/2 



100M Q 
MbhW) 



1/2 



yr(33) 



Therefore, with these units, velocities are given in 
0.6545(M BH (0)) 1 / 2 (lpc/r ) 1 / 2 km/s. 

If the external radius of the disk is fixed at 50 pc, mod- 
els are characterized by three parameters: the masses of the 
seed and of the disk as well as the critical Reynolds num- 
ber. Tabic 1 gives the parameters of the different models 
investigated. The first column identifies the model, the sec- 
ond indicates the disk mass, the third gives the seed mass 
and the fourth gives the critical Reynolds number. 

All computations were performed at the Centre of 
Numerical Computation of the Observatoire de la Cote 
d'Azur (SIGAMM). 



6 



M. Montesinos and J. A. de Freitas Pacheco: The growth of supermassive black holes fed by accretion disks 



Table 1. Disk model parameters 



Model 


Disk mass 


Seed mass 


Reynolds number 




(M Q ) 


(Af ) 




1 


5 x 10 7 


100 


500 


2 


5 x 10 8 


100 


500 


3 


5 x 10 8 


1500 


500 


4 


5 x 10 7 


1500 


500 


5 


5 x 10 7 


1500 


1500 


6 


5 x 10 7 


100 


1500 



1E9 
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Fig. 1. Evolution of the black hole masses. Labels indicate 
the considered models, whose parameters are given in table 
1. 

4. Results 

4.1. The growth of seeds 

The analysis of the numerical solutions obtained for the 
considered models indicate that after ^1.2 Gyr, the disks 
have conserved only ~ 1-2% of their original masses. This 
result depends weakly on the initial mass either of the disk 
or of the seed and are valid for critical Reynolds numbers 
in the range 500 — 1500. About 52-54% of the disk mass 
is accreted by the central BH whereas approximately 45% 
are "lost" as a consequence of the expansion of the outer 
regions, due to the redistribution of angular momentum 
throughout the disk. These values are probably upper limits 
since the star formation process, not included in the present 
approach, may consume a substantial fraction of the gas. 

Figure 1 shows the evolution of the central black hole 
mass for models 1, 2, 4 and 6, since models 3 and 5 have 
a behavior respectively similar to models 2 and 1. A close 
inspection of these curves permits an evaluation of the pa- 
rameter t$o, which measures the timescale required for the 
seed to accrete 50% of the initial disk mass. Table 2 gives 
for each model the mass fraction accreted by the central 
BH after 1.2 Gyr and the parameter in Myr. The mass 
fraction of the disk beyond the initial radius (50 pc) or the 
mass "lost" due to the expansion is also given in table 1 
as well as the parameter t^o, which measures the timescale 
necessary for the disk to lose 40% of its initial mass due 
to the redistribution of angular momentum. Notice that 
for model 6, which has a seed of 100 M Q and a critical 
Reynolds number equal to 1500, the timescale £50 is longer 
than 1.2 Gyr, time interval in which the evolution of the 
disks was followed and, consequently, by that time the BH 
is still growing but at a very slow rate. 



Table 2. Accretion and mass loss parameters: fraction of 
the initial disk mass accreted by the seed after 1.2 Gyr 
(column 2); timescale in Myr for accreting 50% of the initial 
disk mass (column 3); fraction of the initial disk mass "lost" 
by expansion (column 4) and timescale in Myr for "losing" 
40% of the initial disk mass (last column) 



Model Accreted mass tso Mass lost £40 



1 


0.519 


541 


0.466 


157 


2 


0.518 


173 


0.466 


56 


3 


0.518 


173 


0.466 


56 


4 


0.542 


133 


0.442 


94 


5 


0.519 


533 


0.466 


156 


6 


0.468 


1426 


0.462 


385 



Models 1 and 6 differ only by the adopted critical 
Reynolds number. Since the viscous timescale t vis = r 2 /i] <~ 
is related with the accretion timescale, one should ex- 
pect that t 50 for model 6 would be about three times higher 
than for model 1, what is in fact verified from the numbers 
given in table 2. Notice that the timescale £40 for model 6 is 
also larger than that of model 1 by a similar factor. Models 
4 and 5 also differ with respect to the critical Reynolds num- 
ber as models 1 and 6. However, the comparison between 
models 4 and 5 indicates that ratio between the values of t 50 
for these models is slightly higher than the value expected 
simply from the ratio between the critical Reynolds number 
defining each model. The reason is due to fact that these 
models have seed masses considerably higher than models 
1 and 6. A higher seed mass dominates earlier the dynam- 
ics of the inner part of the disk, increasing the accretion 
rate and reducing t$o. Thus, the rapidity at which the cen- 
tral BH accretes 50% of the disk mass depends not only of 
the critical Reynolds number (certainly the main parame- 
ter controlling the accretion process) but also on the initial 
seed mass. Models 1 and 2 have the same seed mass and 
critical Reynolds number but the disk of the latter is ten 
times more massive than that of model 1. More massive 
the disk is, higher the mass density, which leads to a higher 
accretion rate. These effects can be seen in fig. 2, in which 
the evolution of the accretion rate for the different models is 
shown. For all models, in the early evolutionary phase, the 
accretion rate decreases slightly as the inner parts of the 
disk are consumed, except for model 4, in which a small in- 
crease is observed up to t = 8 x 10 yr, followed by a slowly 
reduction of the accretion rate. At the end of this phase, 
the accretion rate decreases abruptly and the growth of the 
BH ceases as well as its activity. 

4.2. The luminosity evolution 

During the early active phase, the BH can be identified 
as an AGN or a quasar. Such a phase, whose timescale is 
given approximately by t 50 lasts for 130-540 Myr, according 
to our models, excepting for case 6, which has a longer du- 
ration due to the smaller accretion rate. These timescales 
are consistent with estimates of the duration of the activity 
based on the statistics of QSOs and AGNs. For instance, 
Marconi et al. (2004) estimated that the activity phase as- 
sociated to SMBH with masses around 4 x 10 8 M Q is about 
450 Myr while the timescale activity related to more mas- 
sive BHs (<~ 10 9 Mq) is considerably shorter, i.e., ~ 150 
Myr, in agreement with our expectations. 
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Fig. 2. Accretion rate evolution for models 1, 2, 4 and 6 



Fig. 3. The bolomctric luminosity evolution for models 1, 
3, 4 and 6 



Figure 3 shows the evolution of the bolometric luminos- 
ity for some of our disk models, which follows, as expected, 
the same trend of the accretion rate history. Maximum 
bolometric luminosities are in the interval 4 x 10 46 — 8 x 10 44 
erg/s, well within the observed range. Notice that a more 
massive disk (~ 4— 5 x 10 9 M©) would produce still higher 
luminosities. The conversion efficiency of gravitational en- 
ergy into radiation derived for our models from the ratio 
e = L/Mc 2 , indicates that such a quantity does not remain 
constant during the disk evolution. In the beginning, it is 
close to the maximum allowed value (e ~ 1/12) and then 
it decreases steadily to values of the order of 10 -3 after 1.2 
Gyr. 

In the early evolutionary phases (t < 100 Myr), the disk 
luminosity exceeds the Eddington value. This question has 
already been extensively discussed in the literature, since 
one may wonder if such a limit should be applied to accre- 
tion disks (see, for instance, Heinzeller & Duschl 2007 for 
a recent discussion on this subject). In fact, the Eddington 
limit expresses the maximum radiative flux able to cross 
the outer layers of a star without destroying its hydro- 
static equilibrium. In general, the gas is supposed to be 
completely ionized and only the scattering of photons by 
free electrons is taken into account in the interaction be- 
tween radiation and matter. Under these conditions, the 
Eddington luminosity depends only on the mass of the star. 
In the case of an accretion disk, the geometry is rather dif- 
ferent and the condition for having equilibrium along the 
z-axis is local. In the region close to the last stable circular 
orbit, tidal forces due to the central black hole balance pres- 
sure gradients along the z-axis of the disk, maintaining the 
hydrostatic equilibrium. If radiation pressure is supposed 
not to surpass gravity, then the following local condition 
should be satisfied 



ad(r) < 



cGM 



BH 

t 



<:f 



(34) 



Notice that the radiative flux is fixed by eq. 19 and, as the 
disk is inflated by the radiation pressure, the energy frac- 
tion advected inwards increases, permitting higher accre- 
tion rates without rising the radiative flux. This "cooling" 
effect is amplified by the photon trapping mechanism, as 
already remarked by Begelman (1978) and Ohsuga et al. 
(2002), contributing to maintain the hydrostatic equilib- 
rium along the vertical axis. 





Model 4 


\ 9.29 Myr 










1 .23 Gyr // 


\ \ 436.6 Myr 





1E12 1E13 1E14 1E15 1E16 1E17 1E18 1E19 1E20 

distance to the black hole (cm) 

Fig. 4. The profile of the aspect ratio H/r for model 4 in 
three different evolutionary stages: 9.29 Myr, 436.5 Myr 
and 1.235 Gyr. 



4.3. Physical properties of the disk 

As mentioned previously, the radiation pressure inflates the 
inner region, producing a slim disk. This effect is illustrated 
in figure 4, in which the profile of the aspect ratio H/r for 
model 4 is shown for three different evolutionary stages. 
At 9.29 Myr after the beginning of the accretion process, 
the radiation pressure is responsible for the increase of the 
scale of height inside a region r < 5 x 10 12 cm. At this 
moment, the effective temperature at the inner radius at- 
tains a value of about 1.8 x 10 6 K. After 436.5 Myr, the 
radiation pressure still affects the inner region, since the ef- 
fective temperature is still of the order of ~ 10 6 K. Near the 
end of the evolution, after 1.23 Gyr, when all the disk mass 
was practically consumed, effects of the radiation pressure 
are no more seen and the scale of height increases with 
the distance as expected. The same behaviour is seen for 
all models having a critical Reynolds number equal to 500. 
For models with higher values of 1Z (models 5 and 6) this 
not happens. Higher Reynolds numbers decrease the radial 
velocity, increasing the local mass density and decreasing 
the parameter j3. Thus, besides the tidal field of the BH, the 
disk self-gravitation also contributes to maintain the hydro- 
static equilibrium in the inner region, avoiding an impor- 
tant inflation of such a zone as in other models with lower 
values of 1Z. 
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Fig. 5. Effective optical depth profiles for model 4 in three 
different evolutionary stages: 9.29 Myr, 436.5 Myr and 
1.235 Gyr. 



The effective optical depth controls the energy evacua- 
tion by radiation along the vertical axis and the local ef- 
fective temperature. In general, the effective optical depth 
increases with the distance, reflecting essentially a decreas- 
ing gas temperature. A maximum is reached and then the 
optical thickness decreases quite rapidly due to the decrease 
of the columnar mass density and the gas recombination. 
This typical behaviour is shown in figure 5 for model 4. We 
can see that as the disk evolves it becomes more optically 
thin in the inner region due to a decreasing columnar mass 
density, consequence of the accretion process. Notice that 
each curve ends at the distance at which the disk becomes 
completely neutral. In this model, the ionized region in- 
creases from a radius of ~ 8 x 10 15 cm up to ~ 2 x 10 17 
cm in the time interval 9.3 — 436 Myr and then decreases 
to ~ 6 x 10 16 cm at the end of our calculations, i.e., 1.24 
Gyr. 

The highest effective temperatures reached in the inner 
region of different models after ~ 9.3 Myr are of the order 
of 8 — 15 millions K obtained respectively for models 1 and 
6. These two models differ only by the critical value of the 
Reynolds number and, as explained above, the former has 
a geometrical thick inner zone while the latter has a thin 
inner zone. Consequently, the advection cooling is more im- 
portant in model 1 than in model 6 and, in spite of having a 
higher viscous dissipation, its central effective temperature 
is lower than that of model 6. The other models after ~ 9.3 
Myr have central effective temperatures of about 2 millions 
K, excepting model 5 which is our "coolest" disk, since its 
effective temperature at that moment is only ~ 400, 000 
K. Models 5 and 6 differ only by the initial BH seed mass 
and have a geometrical thin inner scale of height as already 
mentioned but nevertheless their central effective tempera- 
tures differ by a factor of four. The main reason for such a 
difference is found in their mass profiles. The viscous dis- 
sipation is proportional to the columnar mass density (see 
eq. 20 ) and, at t ~ 9.3 Myr, £ in model 6 is about two or- 



20] ) and, at t 

ders of magnitude higher than in model 5, since the higher 
seed mass of model 5 consumes faster the gas in the inner 
parts of the disk than in model 6. 

Effective temperature maps for models 5 and 6 are 
shown in figure 6. Notice the higher values at the inner 
disk reached in model 6 due to the reason already men- 
tioned. However, the temperatures in the neutral zone are 
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Fig. 6. Effective temperature maps for models 5 (left panel) 
and model 6 (right panel) in three different evolutionary 
stages: 9.29 Myr, 436.5 Myr and 1.235 Gyr. 



quite similar in both models: T ~ 100 K after 9.29 Myr, 
T ~ 500 K after 436.5 Myr and T - 670 K after 1.235 Gyr. 
The temperature in the neutral zones of disk models with 
primordial chemical composition results from the balance 
between viscous dissipation and radiation by molecular hy- 
drogen excited by collisions with H atoms. 

The evolution of the columnar mass density profile de- 
pends essentially on the inward and outward mass fluxes. 
The former is controlled mainly by the critical Reynolds 
number and by the mass of the seed while the latter de- 
pends mainly on the redistribution of the angular momen- 
tum throughout the disk. In other words, in the inner re- 
gions the columnar mass density decreases as the seed grows 
and in the outer regions, the columnar mass density de- 
creases because the disk expands. Consequently, at the end 
of its evolution, the remaining mass of the disk is concen- 
trated in a "torus-like" structure rotating around the cen- 
tral BH . This behaviour is shown in figure 7, in which the 
evolution of the columnar mass density profile for model 
3 is presented. After ~ 1.2 Gyr the "torus-like" structure 
having a diameter of the order of 2 pc is clearly seen. The 
gas temperature in this region is about 2000 K and it ro- 
tates with a Keplerian velocity of about 2640 km/s, con- 
trolled by the central BH mass. It worth mentioning that 
in NGC 1068, the central BH is surrounded by a torus-like 
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Fig. 7. The evolution of the columnar mass density profile: 
time in the different snapshots corresponds to the same 
order as in previous figures. Notice the torus-like structure 
in the last panel, at the end of the disk evolution. 

structure with a diameter of 3.4 pc, detected by the dust in- 
frared emission (Jaffe et al. 2004). The torus-like structure 
is an important element of the so-called "unified model" 
of AGNs and quasars since, according to the direction of 
the line of sight, it determines the visibility (or not) of the 
"broad line region" . It is remarkable that our model is able 
to offer a natural explanation for the formation of such a 
"torus-like" structure as a simple consequence of the disk 
evolution. 

5. Conclusions 

In the present work are reported the results of numerical 
solutions of the hydrodynamical equations describing non 
stationary accretion disks. These simulations aim to inves- 
tigate the fraction of the disk mass transferred to the cen- 
tral BH and to obtain the timescales of such a process. 
Turbulence generated by gravitational instability is sup- 
posed to be the source of viscosity responsible for the an- 
gular momentum transport throughout the disk. The flow, 
self-regulated by such a mechanism, is characterized by a 
critical Reynolds number with values in the range 500-1500. 
The other parameters affecting the accretion process are the 
initial masses of the seed and of the disk itself. 

Observations of bright quasars at z ~ 6 associated to 
SMBH having masses ~ 10 9 M Q imply that their growth by 
accretion processes must occur in timescales less than 900 
Myr. The present investigation suggests that our models are 
able to feed seeds with masses in the range 100 — 1500 M 
in the required timescale. Our numerical simulations in- 
dicate that seeds capture about a half of the original disk 



mass, if the gas consumption by star formation is neglected. 
This conclusion is weakly dependent on the adopted model 
parameters as the initial masses of the BH and of the 
disk or the critical Reynolds number. However, the growth 
timescale, the accretion rate as well as the resulting disk 
luminosity do depend on the adopted value of those pa- 
rameters. Thus, in order to obtain a quasar at z ~ 6 ra- 
diating with a power of ~ 4 x 10 47 erg/s, a disk of few 
10 9 M & is necessary and the accretion process should be- 
gin around z ~ 8. However, it should be emphasized that in 
such massive disks the star formation process is likely to be 
very efficient, which may reduce considerably the amount 
of gas available to feed the black hole. In fact, Shlosman 
& Begelman (1989) concluded from their investigation that 
even a low efficiency in the gas conversion process into stars 
would deplete the disk on a relatively short time scale. 

Although the formation of several circumnuclear disks 
during the history of a galaxy cannot be excluded, it is 
likely that the growth of seeds is a consequence of a unique 
episode. In this case, one should expect that SMBHs grow 
faster than their host halos. In fact, observations suggest 
that the growth of BHs is linked to the evolution of galax- 
ies but not necessarily linked to the growth of their halos. 
Such an "anti-hierarchical" growth is suggested by the fact 
that the co-moving number density of low X-ray luminosity 
AGNs peak at z < 1 while that of high X-ray luminosity 
AGNs peak at z ~ 2 (Ueda et al. 2003; Hasinger, Miyaji 
& Schmidt 2005). Since one should expect that the X-ray 
luminosity is directly related to the accretion rate and, con- 
sequently, to the BH growth, these observations suggest an 
early assembly of the more massive objects. In fact, such an 
interpretation is consistent with our simulations, since mod- 
els with higher seed masses have higher accretion rates and 
higher luminosities. Our computations also show that the 
conversion efficiency of gravitational energy into radiation, 
measured by the ratio L/Mc , does not remain constant 
during the evolution of the disk, varying within the range 

lo^-nr 3 . 

Detection of polarization in the broad emission lines of 
NGC 1068 provided a major argument in favour of the uni- 
fied AGN model (Antonucci 1993). According to this uni- 
fied scenario, AGNs have a central core often displaying 
jets, a inner broad line region (BLR), a narrow line region 
(NLR) and a circumnuclear torus constituted of dust and 
gas. Observation or not of the BLR depends on the inclina- 
tion angle of the line of sight with respect to the plane of 
the torus. Recent observations (Jaffe et al. 2004) indicate 
that the torus is quite compact (few parsecs sized) and hav- 
ing probably a "clumpy" structure. Its origin is still under 
debate (see, for instance, Elitzur & Shlosman 2006) but 
such a structure appears quite naturally in our models as a 
consequence of the disk evolution. The torus-like structure 
results from the matter consumption in the inner regions 
by the central BH and matter "lost" in the outer regions as 
a consequence of the disk expansion by transfer of angular 
momentum. In the transition zone, where the radial veloc- 
ity passes from negative to positive values, the disk material 
remains stagnant and forms the aforementioned structure. 
The dimensions of the resulting torus for model 3, as men- 
tioned above, are quite compatible with the parsec-sized 
torus observed in NGC 1068 (Jaffe et al. 2004). Moreover, 
the Toomre's parameter in the torus satisfies the condition 
Q < 1, indicating that the gas is gravitationally unstable, 
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which could be an explanation for its clumpy structure, 
supporting our formation scenario. 

Our simulations indicate that a substantial fraction of 
the disk remains neutral with the gas temperature in the 
interval 100 — 2000 K most of the time. This late evolu- 
tion of the outer regions of accretion disk could be related, 
for instance, to the molecular "ring" of 2 pc radius ob- 
served around Sgr A* (Gusten et al. 1987). Moreover, the 
physical conditions prevailing in the outskirts of the disk 
are favourable to star formation and could be an explana- 
tion for the presence of massive early-type stars, located in 
two rotating thin disks around the BH situated in our own 
galaxy (Genzel et al. 2003; Paumard et al. 2006). These 
stars have an estimated age of 6 — 8 Myr and the total 
mass under the form of stars and gas in these disks is about 

2 x 10 4 Mq. It is interesting to notice that according to our 
simulations, in order to produce a SMBH of mass of about 

3 x 10 6 Mq, like the one located in the Milky Way centre, a 
disk of initial mass of about 6 x 10 6 Mq is necessary. At the 
end of its evolution, it would remain a mass of about few 
10 4 Mq, consistent with observations. This scenario will 
be investigated in more details in a future paper reporting 
simulations including the possibility of gas conversion into 
stars. 
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